#' Add a 'gages' column to a RouteLink.nc file
#' 
#' @param rlFile Character path+file for the desired Route_Link.nc file
#' @param gageIds Vector of character identifiers for the gages, max length = 15. Must collocate with comIds.
#' @param comIds  Vector of integer comIds for identifying the links (existing in rlFile). Must collocate with comIds.
#' @param newCopyId Character identifier for the new copy of the fullDomFile with specified frxst_pts.
#' @param gageMiss Value which indicates no gage at this link.
#' @param overwrite Logical If the output path/file created from fullDomFile and newCopyId already exists, then overwrite it.
#' 
#' @examples
#' \dontrun{
#' rlFile <- '~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.nc'
#' gageIds <- c('06730200', #BOULDER_CREEK_AT_NORTH_75TH_ST._NEAR_BOULDER_CO
#'              '06730160', #FOURMILE_CANYON_CREEK_NEAR_SUNSHINE_CO
#'              # '06727410', #FOURMILE_CREEK_AT_LOGAN_MILL_ROAD_NEAR_CRISMAN_CO
#'              '06727500'  #FOURMILE_CREEK_AT_ORODELL_CO
#'              )
#' comIds <- c( 42,  #BOULDER_CREEK_AT_NORTH_75TH_ST._NEAR_BOULDER_CO
#'              38,   #FOURMILE_CANYON_CREEK_NEAR_SUNSHINE_CO
#'              # 200,  #FOURMILE_CREEK_AT_LOGAN_MILL_ROAD_NEAR_CRISMAN_CO
#'              82   #FOURMILE_CREEK_AT_ORODELL_CO
#'               )
#' newCopyId <- 'threeRealGagesTEST'
##' AddRouteLinkGage(rlFile, gageIds, comIds, new=identifier)
##' if(FALSE) {
##' library(rwrfhydro)
##' rlFile <- '~/RouteLink_2015_12_15.nc'
##' p <- ncdump("~/nudgingParams.conusPstActive.nc",'stationId',q=TRUE)
##' g <- ncdump(rlFile,'gages',q=TRUE)
##' whGAndP <- which(g %in% p)
##' whGNotP <- which(!(g %in% p))
##' (length(whGAndP)+length(whGNotP))==length(g)
##' g[whGNotP] <- formatC('', width=15)
##' setdiff(g,p)
##' setdiff(p,g)
##' comIds <- ncdump(rlFile,'link', q=TRUE)
##' newCopyId <- 'conusPstActive'
##' gageIds <- g
##'
##' ## NO HI and PR
##' library(rwrfhydro)
##' #rlFile <- '~/RouteLink_2016_02_19_no_HI_PR.nc'
##' rlFile <- '~/WRF_Hydro/TESTING/TEST_FILES/CONUS/WORKING/DOMAIN/RouteLink_2016_02_19_no_HI_PR_goodlakes1260.nc'
##' p <- ncdump("~/nudgingParams.conusPstActive.nc",'stationId',q=TRUE)
##' g <- ncdump(rlFile,'gages',q=TRUE)
##' whGAndP <- which(g %in% p)
##' whGNotP <- which(!(g %in% p))
##' (length(whGAndP)+length(whGNotP))==length(g)
##' g[whGNotP] <- formatC('', width=15)
##' setdiff(g,p)
##' setdiff(p,g)
##' comIds <- ncdump(rlFile,'link', q=TRUE)
##' newCopyId <- 'conusPstActive'
##' gageIds <- g
##' }
#' }
#' @keywords manip IO
#' @concept nudging dataMgmt
#' @family nudging
#' @export
AddRouteLinkGage <- function(rlFile, gageIds, comIds, newCopyId, gageMiss='', overwrite=FALSE) {

  if(length(gageIds)!=length(comIds)) 
    warning("gageIds and comIds do not have same length.", immediate.=TRUE)
  
  if(missing(newCopyId)) warning('A newCopyId required to distinguish the output file with gages.')
  
  rl <- as.data.frame(GetNcdfFile(rlFile, quiet=TRUE))
  gagesInOrig <- any(names(rl)=='gages')
  rl$gages <- formatC(gageMiss, width=15) ## populate with missing

  if(nrow(rl)==length(gageIds)) {
    rl$gages <- gageIds
  } else {
    for(ii in 1:length(gageIds)) 
      rl$gages[which(rl$link==comIds[ii])] = formatC(substr(gageIds[ii], 1, 15),width = 15)
  }

  newDir <- dirname(rlFile)
  origBase <- basename(rlFile)
  newFile <- paste0(newDir, '/', gsub('\\.',paste0('.',paste0(newCopyId,'.')), origBase))
  if(file.exists(newFile) & !overwrite) {
    warning(paste0(newFile, ' exists and overwrite is not specified. Returning.')) 
    return('')
  }
  
  file.copy(rlFile, newFile, overwrite = TRUE)
  
  ncid <- ncdf4::nc_open(newFile, write=TRUE)
  link <- ncid$dim[['linkDim']]
  idDim <- ncdf4::ncdim_def( 'IDLength', '', 1:15)
  gages <- ncdf4::ncvar_def('gages', 'usgsId', list(idDim,link), formatC(gageMiss, width=15))
  if(!gagesInOrig) ncid <- ncdf4::ncvar_add(ncid, gages)
  ret <- ncdf4::ncvar_put( ncid, 'gages', rl$gages)
  ncdf4::nc_close(ncid)
  newFile
}


#' Edit the a frxst pts layer in Fulldom
#' 
#' Note that this function removes any existing frxst_pts by default. The keep keyword keeps existing frxst_pts
#' upto their replacement by ones specifed by gridINds and frxstInds. 
#' 
#' @param fullDomFile Character file/path to the Fulldom file to copy.
#' @param newCopyId Character identifier for the new copy of the fullDomFile with specified frxst_pts.
#' @param gridInds Integer vector of gird indices to assign the integer values in frxstInds.
#' @param frxstInds Integer identifiers for gridInds. 
#' @param keep Logical Keep existing forecast points in the frxst_pts layer. May still be replaced by frxst_pts specfied 
#'              by gridInds and frxstInds.
#' @param overwrite Logical If the output path/file created from fullDomFile and newCopyId already exists, then overwrite it.
#' 
#' @examples 
#' \dontrun{
#' # this example shows how to roughly match frxst_pts to reach gages
#' fullDomFile <- '~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Fulldom_hires_netcdf_file.nc'
#' newCopyId <- 'matchRlGages'
#' gages <- ncdump("~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.threeRealGagesTEST.nc", 
#'                 'gages', quiet=TRUE)
#' link <- ncdump("~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.threeRealGagesTEST.nc", 
#'                'link', quiet=TRUE)
#' whGages <- which(gages != '               ')
#' gageLinks <- link[whGages]
#' names(gageLinks) <- gages[whGages]
#' linkId <- ncdump(fullDomFile, "LINKID", quiet=TRUE)
#' whGagesFull <- plyr::llply(gageLinks, function(gl) which(linkId == gl)[1])
#' EditFrxstPts(fullDomFile, newCopyId, 
#'              gridInds=unlist(whGagesFull),
#'              #frxstInds=as.integer(names(whGagesFull)),  ## fails because these cant be represented as short integers
#'              frxstInds=as.integer(c(500,160,200)), 
#'              overwrite=TRUE)
#' }
#' @keywords manip IO
#' @concept nudging dataMgmt
#' @family nudging
#' @export
EditFrxstPts <- function(fullDomFile, newCopyId, gridInds, frxstInds, keep=FALSE, overwrite=FALSE) {

  newDir <- dirname(fullDomFile)
  origBase <- basename(fullDomFile)
  newFile <- paste0(newDir, '/', gsub('\\.',paste0('.',paste0(newCopyId,'.')), origBase))
  if(file.exists(newFile) & !overwrite) {
    warning(paste0(newFile, ' exists and overwrite is not specified. Returning.')) 
    return('')
  }

  file.copy(fullDomFile, newFile, overwrite = TRUE)
  
  ncid <- ncdf4::nc_open(newFile, write=TRUE)
  frxstPts <- ncdf4::ncvar_get(ncid, 'frxst_pts')
  if(!keep) {
    # currently the specified missing_value is inaccurate, so it's futile to code 
    # up their proper handling
    frxstPts[] <-  as.integer(-9999)
  } 
  frxstPts[gridInds] = as.integer(frxstInds)
  #gages <- ncdf4::ncvar_def('gages', 'usgsId', list(idDim,link), formatC(gageMiss, width=15))
  ret <- ncdf4::ncvar_put( ncid, 'frxst_pts', frxstPts)
  ncdf4::nc_close(ncid)

  newFile
}











